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I show that the basic structure of symplectic integrators is governed by a theorem which states 
precisely, how symplectic integrators with positive coefficients cannot be corrected beyond second 
order. All previous known results can now be derived quantitatively from this theorem. The 
theorem provided sharp bounds on second-order error coefficients explicitly in terms of factorization 
coefficients. By saturating these bounds, one can derive fourth-order algorithms analytically with 
arbitrary numbers of operators. 
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Algorithms for solving diversephysical problems rang- 
ing from celestial mechanics^, El 13 > quantum statis- 
tical mechanics 0, El 0, HI to quantum dynamics [Tol 
EH IT^ . Il3| can all be derived from approximating the 
evolution operator e e ^ T+v ^ in the product form 



n 



e UeT e v^V 



exp 



e(e T T + e v V + ee T v[T, V] 



-e 2 e TT v[T, [T, V}} + e 2 e VT v[V, [T, V}} 



(1) 



with factorization coefficients {ti} and {vi}. Classi- 
cally, every product of the form Q produces a sym- 
plectic integrator for integrating classical equations of 
motion. For solving time-irreversible problems involving 
the diffusion operator, such as the quantum statistical 
trace 0, IE B , H, or the imaginary time Schrodinger 
equation]]! Mill 0, 

one must also insists that 
these coefficients be positive. Since T and V are non- 
commuting operators, the general problem of deriving 
approximations of the form Q beyond second order, re- 
gardless of the sign of the coefficients, is extremely dif- 
ficult. For N > 3, most higher order algorithms can 
only be found by usi ng s ymbolic algebra and numerical 
methods [A ll3L 1 1 8L Il9l |20 | . In this work, I prove a funda- 
mental theorem relating the error coefficients exv, &ttv 
and evTV from which one can deduce previous known 
results quantitatively and construct fourth-order algo- 
rithms analytically with arbitrary N. 

The error coefficients er, ey, erv, s-ttv, and evrv 
in are related to the factorization coefficient {ti} and 
{vi} via 8] 
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r ^2u(si + Si_i)«j , (4) 
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Satisfying the primary constraints ej- = 1 and ey = 1 
implies that sn = 1 and 



U\ = 1. 



(7) 



For {ti} > 0, JSJ is a quadratic form in Ui whose min- 
imum can be determined subject to linear constraints 
(0 and Q. This then leads to an inequality relating 
e-TV, &ttv 1 and evrv sufficient to prove that the general 
product iQJ cannot be corrected (to be explained below) 
beyond second order with positive coefficients^- How- 
ever, the exact minimum was not determined because the 
required sum, first appeared in Suzuki's workpl). 
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Sg), 
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was not known in closed form. The inequality was there- 
fore weak, excluding the possibility of being equal. Sur- 
prisingly, a closed form exists and was found in Ref . p^ ] , 



i=l 



(9) 



The minimum was then determined, but without explic- 
itly incorporating the constraint Q. Recently, it was 
realized [23j that constraint (J7J) can be enforced without 
affecting any of the equations J3J, and JSJ by sim- 
ply setting t\ = 0. The resulting minimum is then only 
true for algorithms whose first operator is e VieV . Since 
classically this operator updates the velocity (momen- 
tum) variable, the constraint J7J) dictates that the min- 
imum derived in Ref. |22j| only holds for velocity-type al- 
gorithms. By interchanging T <-> V and {ti} «-> {vi} in 
all of the above, the constraint ey = 1 now dictates that 
Vi = and another minimum holds for position type al- 
gorithm whose first operator is e tieT . One is finally able 
to state the exact relationship between eyy, &TTVi an d 
evrv directly in terms of either {ti} or {vi} correspond- 
ing to either velocity or position-type algorithms. By 
constructing integrator whose error coefficients are pre- 
cisely at the quadratic minimum, the condition for being 
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fourth-order can be directly stated, and easily solved for, 
in terms of {ti} or {vi}. One is then able to construction 
fourth-order integrators analytically for arbitrary N as 
it was done in Ref..23j. The current theorem provided 
sound theoretical support and unified derivation of re- 
sults obtained in Ref . I2H . 



II. THE THEOREM 

The constrained minimum of the quadratic form in (J5J 
can be obtained by the method of Lagrange multiplier. 
Since this has been worked out in details in Ref.|8j] (but 
for a much weaker goal), we will just summarize the re- 
sults. For t\ = and {ii>i} > 0, minimize 

1 N ( N \ ( N 

F = 2 E ^ 2 - Ai E -ME 

i=l \i=l / Vi=l 



with respect to Ui gives, 

Ui — Ai + A 2 (s.i + Si-i). 




(11) 



Substituting this back to satisfy constraints @ and Q 
determines Ai and A 2 : 



where Sg is given by ©. Since < Sg < 1 for {ti>i} > 
and ex = 1, the second form shows that the LHS of 
is strictly negative. Note that t\ = does not prevent the 
algorithm from being completely general. Nothing stops 
us from considerating algorithms with v\ = 0, in which 
case, the result will be a position-type algorithm. This 
part of the theorem simply regard {ti} as independent 
variables. 

By interchanging T ^ V and {ti} <-> {v^ in QJ, 
the error coefficients changes respectively, erv — * —£tv, 
&ttv —zvtv an d evrv — * —&ttv- Making the sub- 
stitution in (|15|l yields, 

Theorem, Part B: For v\ = 0, and {vi>i} > 0, the 
error coefficients for the product of operators 
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obey the inequality, 
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(17) 
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Ai + A 2 — — + exv, 



(12) 



Ai + A 2 + 3A2 = - + e T v + 2gttv, (13) 
o 



where J3i)Li ^i( s i + 5 i-i) 2 =1 + 3- The minimum of the 



quadratic form is then 
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i[(A 1+ A 2 ) 2 + ff A2], 
i(i + e Ty ) 2 + i-(2e T TV-i) 2 . (14) 
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Setting the LHS of JHJ greater or equal (this is the most 
important point, the main contribution of this work) to 
F mm gives, 

Theorem, Part A: For t\ = and {i^>i } > 0, the error 
coefficients for the product of operators in obey the 
inequality, 



e VTV < 
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24 2 1 V 1 - <Sg v 
or, after slight arrangement, 
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evTV + 2 6tv ~ eTTV 



I ' / 1 - J • (I 5 ) 



" -24^"T^ (eTTy -^ 3)2 ' (16) 



or in the form 
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~ 24 y 
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where the corresponding Sg' is given by 



(20) 



Again H19|l shows that the LHS is strictly positive. We 
will regard Hit) [I and (|19fl as fundamental statements of 
our theorem. To explain this, we need to mention sym- 
plectic corrector (or process) algorithms |3 EE E3- 

If p denotes an approximation to e e ( T+v ) of the prod- 
uct form m, then p is "correctable" if 



SpS~ 



(21) 



is correct to higher-order in e for some operator S also of 
the general form but with no sign restriction on its 
factorization coefficients pil I2M 126^ . If p is correctable, 
then its trace, equal to the trace of p', will be correct 
to higher order in e. This is important for calculating 
the quantum statistical trace of an approximate density 
matrix, as in path integral Monte Carlo calculations^ El 
II El The criterion for p to be correctable to at least 
third-order in e isla 



1 1 2 n 
ZVTV + 2 e TV ~ e TTV — u 



(22) 
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However, if {ti} > 0, then (|16|l shows that this is not 
possible. And if {wj} > 0, then (|19fl also shows that 
this is not possible. Our theorem states precisely, how 
forward symplectic integrator of the product form, con- 
sisting of only operators T and V , cannot be corrected 
beyond second order. 

A much weaker form of this theorem, that the LHS 
of (|22l) cannot be zero, has been proved previously by 
Chiny|, and by Blanes and Casas|23 using a very dif- 
ferent method. The current theorem is much sharper, 
stating the precise amount by which the correctability 
condition i)22[l is being missed, when {ti} > 0, and when 
{Vi} > 0. 

Two main corollaries: 1) It is easy to force erv = 0; 
all odd-order error terms will vanish if we simply choose 
factorization coefficients that are left-right symmetric in 
(JTJ or i|17|) . If ztv — 0, then the correctability criterion 
is just errv — eyrv- However, (|Pd|> and (|T§)) both show 
that there is an unbridgeable gap between the two co- 
efficients; they can never be equal. In particular, they 
can never both e qual to zero. This corollary is the Sheng- 
Suzuki Theorem |2ll l28j : there cannot be factorization 
algorithms of the form Q with positive coefficients be- 
yond second order. Again our current theorem is more 
quantitative in showing that if {ti} > 0, then the gap is 
given by H16|) and if if {vi} > 0, then the gap is given 
by |Q . 2) If both erv and exrv are zero, then l(Po|) 
implies that 



evrv < 



1 Sg 
'241-<V 



(23) 



and can only vanish if Sg — 0, requiring at least one U to 
be negative. If both ctv and evTV are zero, then 1|19|) 
implies that 



1 



Sg' 



24 1 - 



(24) 



and can only vanish if Sg' = 0, requiring at least one Vi 
to be negative. This corollary is the Goodman-Kaper 
theorem[2£j: beyond second order, factorization algo- 
rithms of with only operators T and V must have at 
least a pair of negative coefficients (tfe, Vk). Our current 
theorem is again much more quantitative with symmetric 
forms lH3) and (|2ljl . 



III. CONSTRUCTING FOURTH-ORDER 
ALGORITHMS 

Since all odd-order error terms vanish with left-right 
symmetric coefficients, fourth-order algorithms can be 
obtained by forcing both errv and evrv to zero. Let's 
consider first velocity-type algorithms described by Part 
A of the theorem. When e-rv and exrv are both zero, 
the bound for eyrv <|23f) is the actual error coefficient for 
algorithms with Ui given by pif). corresponding to 



-^2(U + ti+i), 



(25) 



Vi = ~ + A 2 (l - t 2 ) and v N = ^ + A 2 (l - t N ), (26) 



with A 2 given by (H2) and fHjl . 



A? 



1 1 



21- Sg' 



(27) 



Ea. (|25() is true for all algorithms whose quadratic form 
is stationary with respect to Ui. Thus the equal sign in 
(12311 holds even for negative ti. A fourth-order algorithm 
results if we choose a left-right symmetric set of {ti} with 
t\ = such that ey = 1 and Sg — 0. For example, for 
N = 6, we can choose t 6 = i 2 , t 5 = < 3 . The constraints 



2/ 2 



2t 3 + U = 1 



24 + 2ti + 4 = 

can be solved by setting i 2 = 0^3, giving 



t 4 - -2 1 / 3 (1 



t 3 , 



t* = 



2(l + a)-2 1 /3(i + a 3)i/3- 



(28) 



(29) 



(30) 



The case of a = reduces back to the well known Forest- 
Ruth integrator!^. For Sg = 0, coefficients Vi given by 
(|2SJl-(|2I3) are linearly related to {ti}. For position-type 
algorithms, one can simply exchange operators T «-> V 
and their coefficients {ti} <-> {v^. For further examples 
of constructing this type of algorithms, see Ref . [23| . 

Instead of forcing evrv to vanish on the RHS of ^ . 
one can simply move the entire operator 

eMe 3 e V Tv[V,[T,V}}) 

back to the LHS, and combine the commutator 
£ c-vtvWi [T, V]} symmetrically with one or more oper- 
ator eviV. For T =_p 2 /2, both classically 13 1| and quan- 
tum mechanically |32|, [V, [T, V]] corresponds to an addi- 
tional gradient force or potential similar to V. By doing 
this, one gets around the Sheng-Suzuki theorem in pro- 
ducing fourth-order forward (ti>i > 0) symplectic inte- 
grator by factorizing e e ( T+v ) through an additional op- 
erator [V, [T,V]]. This results in a far richer family of 
algorithms since any set of symmetric coefficient {ii>i} 
(regardless of sign) satisfying ex = 1 will now yield a 
fourth-order algorithm. For example, for N = 4, taking 
t% = t 3 = ti = 1/3 produces 



1 



Vi = Va 



V2 = V 3 



evrv 



1 

192' 



(31) 



which is forward algorithm 4D[l2|- But one can also take 
t 2 =U = 2/3, t 3 = -1/3, giving 



1 



I'l = l'4 



V 2 = v 3 



evTV 



5 

96' 



(32) 



This also illustrates that the Goodman-Kaper theorem 
no longer holds if one includes [V, [T, V}] in the factor- 
ization process. More examples of deriving velocity- type 
gradient algorithms are given in Ref. |23l| . 
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To construct position-type algorithms, we invoke Part 
B of the theorem. Since we have a preference for keeping 
the commutator [V, [T, V}], we must set erv = &ttv = 
in ljT%|) and solve for eyrv, 



evrv 



12 



(33) 



where we have picked the solution which vanishes with 
Sg'. Under the interchange T «-> V and {ij} <-> {vi}, the 
corresponding equation l|13|) for A2 reads 



1 ^ >\ 1 9 

- + g A 2 = - - 2ey T y, 



from which one deduces 

A * 



1 



(34) 



(35) 



Thus any set of symmetric {v^} with v± = Q and ey = 1 
will produce a fourth-order algorithm via 



U = -\z(vi + v i+ i), 



(36) 



with A2 and eyrv given by i|35|l and (|33|l . For iV 

«2 = ^3 = v a — 1/3, this gives 



h —ti — 



*2 = <3 = 



1 



2V2' 



and 



s-vtv 



:V2 



(38) 



(39) 



More examples of deriving position-type forward algo- 
rithms can be found in Ref.[23j. 

In conclusion, I have presented a fundamental theorem 
on symplectic integrators from which fourth-order algo- 
rithms of arbitrary length can be constructed for solving 
diverse physical problems. Other important properties 
of symplectic integrators can also be deduced from this 
theorem. 

I thank Drs. Blanes and Casas for a comment on 
Ref.[22| which triggered my understanding on the sub- 
ject. This work is supported, in part, by a National Sci- 
ence Foundation grant, No. DMS-0310580. 



h = - + A 2 (l - v 2 ) and t N = ~ + A 2 (l - v N ), (37) 
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